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We propose an experiment to obtain the phase diagram of the fermionic Hubbard model, for any 
dimensionahty, using cold atoms in optical lattices. It is based on measuring the total energy for a 
sequence of trap profiles. It combines finite-size scaling with an additional 'finite-curvature scaling' 
necessary to reach the homogeneous limit. We illustrate its viability in the ID case, simulating 
experimental data in the Bethe-Ansatz local density approximation. Including experimental errors, 
the filling corresponding to the Mott transition can be determined with better than 3% accuracy. 
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Bose-Einstein condensation of trapped atoms [l[ has 
led to rapid growth of research into such systems. One 
exciting suggestion Q is that they could be used to simu- 
late condensed matter, with the atoms playing the role of 
conduction electrons, and a periodic potential provided 
by a laser standing wave, or 'optical lattice'. The advan- 
tages of such a scheme are manifold: the system is highly 
tunable, including the effective strength of atom-atom 
interactions Q; there is no naturally occurring disorder 
in the lattice; and quantities such as the total energy, 
which cannot normally be measured in condensed mat- 
ter systems, are accessible Several familiar phases of 
matter have already been seen: the Fermi liquid [5| , the 
Mott insulator 0], vortex matter 0], and BCS conden- 
sates 0. So far, however, most observations have been 
(a) rather qualitative in nature, and (b) of phases that 
were already well understood. In this Letter, we propose 
a scheme to answer quantitatively an outstanding open 
question of condensed matter physics: what is the phase 
diagram of the fermionic Hubbard model in two or three 
dimensions? 

To answer this question we need to address the fact 
that atom traps are inhomogeneous. Because of that, 
phase transitions are 'blurred' 0]. This is not simply an 
issue of finite size: the system is not in a thermodynamic 
phase at all, even in the thermodynamic limit (see be- 
low). Thus, the notion of a phase diagram for the atomic 
gas in the optical lattice is not rigorous. Nevertheless, 
by focusing on the total energy and carrying out a se- 
quence of measurements for traps of different sizes and 
shapes the phase diagram of the model of interest may 
be deduced, as we shall show. 

To achieve this, we need to link quantitatively the en- 
ergy of the model system (atoms in an optical lattice) to 
that of an infinite d-dimensional Hubbard lattice. One 



obvious issue is size: condensed matter systems typically 
have N ~ 10^^ conduction electrons, whereas N ~ 10^ 
is more typical in the optical lattice case. In computer 
simulations, this is dealt with by finite-size scaling [lo| : 
systems of various sizes are simulated, and the results 
extrapolated to the thermodynamic limit. A similar pro- 
cedure can be employed in optical lattice experiments, as 
we describe below. 

Another difference between solid-state and cold-atom 
systems is the nature of the confining potential. In the 
former, it is usually a hard-wall box; in the latter, it is 
a power-law (usually harmonic) potential. The thermo- 



dynamic limit mentioned above is singular [111 . |12| | : even 
though it corresponds to vanishing trapping potential, 
the results obtained in that limit depend on the shape of 
the trap. Thus finite-size scaling alone docs not yield the 
desired bulk energy, and a second kind of scaling is called 
for: an extrapolation from the properties of the power- 
law-trapped system to those of the hard-wall-confined 
one. We call this finite- curvature scaling. 

Thus we consider the family of potentials given by 
[l^ V{x.) = tJ2i=i Ixi/Lil", where Li are characteris- 
tic length-scales of the trap, and t is a reference energy 
scale (see below), a — 2 corresponds to harmonic trap- 
ping, while the a —>■ oo limit creates a hard-wall box of 
volume 2'^ Hi Li. The d = 1 version of this potential is 
shown in Fig. [Ha). 

In practice, the achievable values of a are likely to be 
the even integers, because the true confining potential 
will generically be harmonic x^) near its minimum. 
To change this one will have to superpose two laser beams 
to cancel the harmonic terms; for Gaussian beams, the 
next term will be quartic ('^ x^). Additional laser beams 
could be used to cancel both the x"^ and x^ terms, leaving 
an x^ potential, and so on. This demands controlling the 



FIG. 1: (a) The trap potential V{x) for d — 1 and various 
values of a. (b) Dependence of total energy on system size for 
a d = 1 optical lattice loaded with non-interacting, spin- 1/2 
fermions. The filling is fixed at f — 1. The trap exponents 
are a = 2 (4-), 4 (*), 6 (□), 8 (■) and oo (x). Inset: the 
same quantity in the thermodynamic limit plotted against 
1/a. The dashed-dotted lines mark the thermodynamic limit 
result for a = oo: Ea/2Lt = -(4/7r) sin(7r//2). 

width of each Gaussian beam independently of its wave- 
length, which is within current experimental capabilities 
|13l . Il4| . Indeed, the realisation of a confining potential 
with a term V{x) ^ has already been reported [15]. 

The experiments we propose are realisations of the fol- 
lowing inhomogeneous Hubbard model 0, [l^ '■ 

H = -tJ2 (cLcjo + H.c.) +[/^c]^CjtsVji 

+ ^F(x,)4c,„ (1) 

jo- 

where cj^ creates a fermionic atom with spin a on site 
j. The terms in ([1]) represent respectively the kinetic 
energy, the potential energy of atom-atom interactions, 
and the potential energy due to the trap. The usual 
Hubbard model, whose phase diagram we wish to deduce, 
is recovered when V{yij) = everywhere. The atoms 
must have two 'spin' states (t and |); in practice these 
will be two hyperfine eigenstates. 

The model ([1} has four dimensionless parameters: U jt, 
the strength of the on-site repulsion; a, the trap ex- 
ponent; A/", the number of sites; and /, the filling. 
These last two are somewhat delicate, since in the pres- 
ence of a power-law trap the effective number of sites 
in the lattice becomes energy-dependent. We thus de- 
fine A/" as the number of sites in the a ^ oo limit, i.e. 
TV = Y^^^{2Li / a) , and / as N /N. The thermodynamic 
limit is 

N oo, N ^ oo, with / = constant (2) 

for given U/t and a, coinciding with the usual definition 
for a — > oo and with Ref. [l3| for a = 2. In the limit 
the trapping potential term in ^ vanishes identically. 



FIG. 2: (a) Same as Fig. [TJb), but with interaction strength 
U/t = 4, and with L/a = 400 data used for the inset. In 
this case the thermodynamic limit result for a = oo is _E = 
—0.57372937. (b) Solid data points and line: the q = 6 data 
from Figs. [TJb) (lower curve) and[2ja) (upper curve). Open 
data points: the same, but with random errors of ±5% added 
to the values of E and /. Dashed line: least-squares fit to 
these 'noisy' data. The resulting error in the thermodynamic 
limit of the energy is < 0.6%. 

To make the case for finite-curvature in addition to 
finite-size scaling, we consider first the non-interacting 
situation. For large enough Af, the density of states 
(DOS) obeys the scahng law p (e) = (f ) . The deriva- 
tion of G{x), which depends implicitly on d and a, has 
been outlined in Refs. [Hill^l and its behaviour for x ^ 1 
has been verified experimentally in the ID case for a — 2 
[H. If we define (f) = J!^^* dxG{x)x, and fix the 

Lagrange multiplier p hy f — j^!^ dxG{x), the ground- 
state energy obeys the scaling law 

E = tMF^ if) . (3) 

Thus Fa{f) is the energy per site in the thermodynamic 
limit ^ in units of t. This allows its determination from 
results for finite size systems by extrapolation. 

Fig. [Ijb) shows some d = 1 total energy curves. The 
energy per site has converged essentially perfectly for sys- 
tem sizes M ^ 200 or larger. This makes it possible to 
obtain Fa(x) from experiments on reasonably small sys- 
tems. However, the scaling function Fa(x) itself is not 
universal, but depends on the trap exponent, a (see fig- 
ure). An additional extrapolation in a is therefore nec- 
essary to obtain the a — ^ oo behaviour (see inset). 

We now turn to the case of interest, U/t > 0. For 
a — > oo and in the thermodynamic limit the model 
^ is known to exhibit at least two distinct phases: a 
metallic one at low or high filling, and an insulating one 
at (and perhaps also around) f — I- The ID Hubbard 
model has been solved exactly [l3|: it is metallic except 
at precisely / = 1, where it is insulating for all J7 > 0. 
But the behaviour of the d > 1 models remains uncertain; 
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for example, it is not clear whether the insulating phase 
occupies just the line / = 1, or some finite region of the 
phase diagram. Moreover, there may be additional com- 
peting phases, including d-wave superconductivity, itin- 
erant antiferromagnetism, and phases with a distorted 
Fermi surface (20| . 

Our proposal for investigating this experimentally 
is as follows: (a) Obtain experimental values for the 
total energy for several values of the filling fraction, 
trap power-law, system size, and interatomic interaction 
strength. This builds up data points from the func- 
tion E{f,U,t,a,N'). (b) Obtain the values of E in the 
thermodynamic limit ([2]) by extrapolation, aided by the 
finite-size scaling law ([3]), where Fa{f) now additionally 
depends on U/t. (c) Numerically differentiate the curves 
thus obtained with respect to the filling /, to obtain the 
chemical potential n{f,U,t,a) = tF^{f,U/t). (d) Take 
the a — » cx) limit by a second extrapolation. This is the 
finite-curvature scaling, and results in a function /i/t of 
/ and U/t. (e) By finding the lines in the f — U/t plane 
where fi/tis discontinuous, construct the phase diagram. 

To demonstrate our method in the one-dimensional 
case, where the answer is already known [l9j, we have 
used density- functional theory within the Bethe-Ansatz 
local-density approximation (BA-LDA) '2l'| to simulate 
the experimental data for step (a). The BA-LDA is 
known to predict ground-state energies of spatially in- 
homogeneous systems with an accuracy of a few percent. 
Furthermore, its accuracy increases as the length-scale of 
the inhomogeneity increases, making it particularly suit- 
able to study the approach to the thermodynamic limit 

Examples of the resulting total energy curves are 
shown in Fig. HJa) for U/t — 4. Note that a scaling 
law of the form ([3]) is still obeyed; indeed, Q can be 
obtained using only the requirement that, for large Af, E 
is extensive. We find that to achieve good convergence 
M ~ 200 sites suffice for a < 8. Experimentally, this 
would require trap frequencies ^ 10 Hz in the ^Li case, 
which is in the currently achievable range. Similarly fast 
convergence is found for larger U/t = 6, 8 (not shown). 

This rapid convergence has the important effect of ren- 
dering the method rather insensitive to experimental er- 
ror. For example, we show in Fig. E^b) the a = 6 curves 
from Figs.mb) and[2Ja), but with random errors of ±5% 
added to the values of / and E. Though these data 
look very noisy, all points for L/a > 30 are essentially 
random scatter around the same value; thus the central 
limit theorem applies, and the overall error may be re- 
duced by increasing the number of measurements. The 
dashed lines are least-squares fits to the noisy data for 
30 ^ L/a ^ 100; their differences from the noise-free val- 
ues are 0.17% (non-interacting case) and 0.61% (interact- 
ing case). Thus a serious error (cumulatively ~ ±8% or 
worse, as Fig. [IJb) shows) in the individual data points 
becomes a negligible one (< 0.6%) in the resulting value 




FIG. 3: Energy per site (a) and numerical estimate of the 
chemical potential (b) vs / for U/t = 4, L/a = 400, a = 2, 4, 
6, 8 and co (lines with no symbols). In the inset in panel (b), 
the filling fmax where fi{f) has maximum slope is plotted as 
a function of (circles). The parabola interpolating the 
points coming from a = 4, 6 and 8 (dashed line) extrapolates 
to fmax{0) = 1.025, indicating that even those relatively small 
trap exponents provide a good estimate of the critical filling. 

for the thermodynamic limit of the energy. 

In Fig. [3l^a) , we plot energy per site as a function of / 
for a = 2, 4, 6, 8, oo and system size Af — 800 (effectively 
in the thermodynamic limit). This must be intepreted 
as a plot of the extrapolated values of E vs nominal val- 
ues of /. The error discussed above (coming from the 
uncertainty in measuring each individual value of E, on 
the one hand, and the fact that the actual filling differs 
from the nominal one for each particular measurement, 
on the other) is smaller than the size of the symbols. The 
cusp indicating the phase transition becomes sharper as 
a is increased, and tends towards its hard-wall position 
at / = 1 as a — > oo, but for finite a it is not a singular- 
ity. Numerically differentiating these curves with respect 
to /, shown in Fig. [3Kb), makes the cusp more apparent 
and shows that there is a single phase transition at / sa 1. 
The inset demonstrates a simple method that allows us 
to obtain the critical value of / within 3% (or ~ 7% if 
only the data points with a < 6 are used). Clearly, the 
same method can also be applied in d = 2 and d = 3, 
where the exact answer is not known. 

The implementation of our proposal requires a signif- 
icant degree of automation. However, preparation-to- 
measurement times of less than a minute are already 
routinely achieved, and complete automation is being 
developed for drop-tower-based zero-gravity experiments 
[2^ . so this requirement is not unfeasible. One could 
alternatively avoid the finite-curvature issues by creat- 
ing a 'flat' potential, using the endcap-beam technique 
2J]. However, Fig.[2Ka) shows that the convergence with 
system-size is much slower in this case, which would lead 
to appreciably larger errors. 

In conclusion, we have presented an experimental 
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scheme to determine the Hubbard model phase diagram 
by measuring the total energy of optical lattice systems. 
It requires a non-trivial curvature extrapolation in addi- 
tion to the usual finite-size one. We have demonstrated 
its validity in the d = 1 case by an approximate numerical 
many-body calculation to simulate experimental data, in- 
cluding an analysis of the effects of realistic experimental 
errors. We conclude that the fermionic Hubbard model's 
phase diagram can be determined quantitatively using 
the experimental procedure that we propose. 

We stress that the experimental use of our method 
would yield the location of the phase boundaries of the 
Hubbard model as functions of interaction strength, [/, 
and filling, /. This would be a significant step forward 
in d — 2 and d = 3, where standard computational ap- 
proaches encounter serious difficulties and relatively little 
is known for certain about the phase diagram. It would 
not, however, permit the identification of the phases; for 
this, further experiments would be required. An obvi- 
ous choice would be to probe correlation functions, e.g. 
via the dynamical structure factor S{a^) using stim- 
ulated emission in a two-laser set-up [251 ) or shot-noise 
techniques [1^. Alternatively, one could exploit the 
fluctuation-dissipation theorem to derive them from the 
total energy; we will discuss this elsewhere. 
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